# TODO: Add comment
# 
# Author: broken
###############################################################################
library(raster)
library(base)
library(gstat)
library(RPostgreSQL)
dr=dbDriver("PostgreSQL")

inputPath= "E:/FIMO/APOM/data/pm/2014/"
resultPath="E:/FIMO/APOM/data/pm/result/2014/"
outputPath="E:/FIMO/APOM/data/pm/newvalidate/2014/"

listFile = list.files(path = inputPath,pattern=".tif$",full.names = FALSE,recursive=TRUE)

listFile=substr(listFile,1,30)

listSize=length(listFile)
listTime=c(1:listSize)

sCor=c(1:listSize)
sRMSE=c(1:listSize)
sRE=c(1:listSize)

oCor=c(1:listSize)
oRMSE=c(1:listSize)
oRE=c(1:listSize)

uCor=c(1:listSize)
uRMSE=c(1:listSize)
uRE=c(1:listSize)

numberStaion=c(1:listSize)

con=dbConnect(dr,dbname="fimo",host="192.168.3.4",port=5432,user="rasdaman",password="rasdaman")
timeQuery="SELECT aqstime FROM apom.satresampmod04 where filename like '%"
selectStr="SELECT stationid,name,st_asText(location),avg(pm25) FROM apom.grdpmhiscem_data,public.grdstation WHERE grdpmhiscem_data.stationid = grdstation.id AND grdstation.id IN (1,2,4,5) AND  aqstime between '"
andStr="' and '"
groupStr="'  group by stationid,name,location having avg(pm25)>0"


for(i in 1:listSize){
	
	skFile=paste(resultPath,listFile[i],"_sk.tif",sep="")
	okFile=paste(resultPath,listFile[i],"_ok.tif",sep="")
	ukFile=paste(resultPath,listFile[i],"_uk.tif",sep="")
	
	ogRaster=raster(paste(inputPath,listFile[i],".tif",sep=""))
	skRaster=raster(skFile)
	okRaster=raster(okFile)
	ukRaster=raster(ukFile)
	
	qr=paste(timeQuery,listFile[i],"%'",sep="")
	resultSet=dbSendQuery(con,qr)
	record=fetch(resultSet,n=1)
	asqTime=substr(record[1,1],1,30)
	myTime=strptime(asqTime,format="%Y-%m-%d %H:%M:%S")
	startTime=myTime-1800
	endTime=myTime+1800
	
	selectQuery=paste(selectStr,startTime,andStr,endTime,groupStr,sep="")
	pmRs=dbSendQuery(con,selectQuery)
	pmRecord=fetch(pmRs,n=-1)
	numberRow=nrow(pmRecord)
	numberStaion[i]=numberRow
	if(numberRow>0){
		ogValueArray=c(1:numberRow)
		skValueArray=c(1:numberRow)
		okValueArray=c(1:numberRow)
		ukValueArray=c(1:numberRow)
		
		for(j in 1:numberRow){
			locationStr=pmRecord[j,3]
			xLocation=substr(locationStr,7,13)
			yLocation=substr(locationStr,15,nchar(locationStr)-1)
			
			ogValue=extract(ogRaster,cbind(as.numeric(yLocation),as.numeric(xLocation)))
			skValue=extract(skRaster,cbind(as.numeric(yLocation),as.numeric(xLocation)),buffer=25000,fun=mean)
			okValue=extract(okRaster,cbind(as.numeric(yLocation),as.numeric(xLocation)),buffer=25000,fun=mean)
			ukValue=extract(ukRaster,cbind(as.numeric(yLocation),as.numeric(xLocation)),buffer=25000,fun=mean)
			
			ogValueArray[j]=ogValue
			skValueArray[j]=skValue
			okValueArray[j]=okValue
			ukValueArray[j]=ukValue
			
		}
		pmRecord$ogColumn=ogValueArray
		pmRecord$skColumn=skValueArray
		pmRecord$okColumn=okValueArray
		pmRecord$ukColumn=ukValueArray
		outFile=paste(outputPath,listFile[i],"_validate.csv",sep="")
		write.csv(pmRecord,file=outFile)
		
		#caculate correlation
		sCor[i]=cor(pmRecord$avg,pmRecord$skColumn)
		sRMSE[i]=sqrt(sum((pmRecord$avg-pmRecord$skColumn)^2)/nrow(pmRecord))
		sRE[i]=sum(abs(pmRecord$avg-pmRecord$skColumn)/pmRecord$avg)/nrow(pmRecord)
		
		oCor[i]=cor(pmRecord$avg,pmRecord$okColumn)
		oRMSE[i]=sqrt(sum((pmRecord$avg-pmRecord$okColumn)^2)/nrow(pmRecord))
		oRE[i]=sum(abs(pmRecord$avg-pmRecord$okColumn)/pmRecord$avg)/nrow(pmRecord)
		
		uCor[i]=cor(pmRecord$avg,pmRecord$ukColumn)
		uRMSE[i]=sqrt(sum((pmRecord$avg-pmRecord$ukColumn)^2)/nrow(pmRecord))
		uRE[i]=sum(abs(pmRecord$avg-pmRecord$ukColumn)/pmRecord$avg)/nrow(pmRecord)
	
	
		#print(paste("Complete ",listFile[i],sep=" "))
		
	}else{
		#print(paste(listFile[i],"have no data ",sep=" "))
	}
	
	
}
statis_result=data.frame(listFile,numberStaion,sRMSE,sRE,oRMSE,oRE,uRMSE,uRE)
outFile=paste(outputPath,"_result_2014.csv",sep="")
write.csv(statis_result,file=outFile)
print("Finish All")



